/*****************************************************************************
 * FILE: train-main.c
 * AUTHOR: William Stafford Noble
 * CREATE DATE: 2/11/99
 * PROJECT: SVM
 * COPYRIGHT: 1999-2001, Columbia University (see ../doc/copyright.html)
 * VERSION: $Revision: 1.53 $
 * DESCRIPTION: Main routine for compute-weights and related programs.
 *****************************************************************************/
#include "kernel.h"
#include "rdb-matrix.h"
#include "compute-weights.h"
#include "classify.h"
#include "kfd.h"
#include "kernel-pca.h"
#include "compute-kernel.h"
#include "class-array.h"
#include "linear-algebra.h"
#include "matrix.h"
#include "array.h"
#include "utils.h"
#include <time.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>

#ifndef KERNEL_PCA
#define KERNEL_PCA 0
#endif

#ifndef KFD
#define KFD 0
#endif

#ifndef COMPUTE_WEIGHTS
#define COMPUTE_WEIGHTS 0
#endif

/* Columns in compute-weights output. */
#define TRUE_CLASS_COLUMN 0
#define WEIGHT_COLUMN 1
#define TRAIN_CLASS_COLUMN 2
#define TRAIN_DISC_COLUMN 3
#define HOLD_OUT_CLASS_COLUMN 4
#define HOLD_OUT_DISC_COLUMN 5
#define NUM_COLUMNS 6

/*****************************************************************************
 * Initialize a set of SVM weights, either by reading from a given
 * file or setting them all to INITIAL_WEIGHT.
 *****************************************************************************/
#define INITIAL_WEIGHT 0.5

ARRAY_T* initialize_weights
  (BOOLEAN_T format_line,
   int       num_examples,
   char*     initial_filename)
{
  FILE*         initial_file;
  RDB_MATRIX_T* initial_weights_rdb;
  ARRAY_T*      weights;
  int           i_weight;

  /* Initialize the weights. */
  if (initial_filename != NULL) {

    /* Open the file. */
    if (open_file(initial_filename, "r", TRUE, "initial weights",
		  "initial weights", &initial_file) == 0) {
      exit(1);
    }

    /* Read the initial weights into an RDB matrix. */
    initial_weights_rdb = read_rdb_matrix(format_line, NULL, initial_file);

    /* Extract the weights from the RDB matrix. */
    weights = get_matrix_column(WEIGHT_COLUMN, 
				get_raw_matrix(initial_weights_rdb));

    /* Make all the weights positive. */
    for (i_weight = 0; i_weight < num_examples; i_weight++) {
      set_array_item(i_weight, abs(get_array_item(i_weight, weights)),
		     weights);
    }

    free_rdb_matrix(initial_weights_rdb);

  } else {

    /* All the weights are set to a single value. */
    weights = allocate_array(num_examples);
    init_array(INITIAL_WEIGHT, weights);
  }

  return(weights);
}


/****************************************************************************
 * Put an RDB header on the output file.
 ****************************************************************************/
static void print_weight_header
  (KERNEL_T* kernel,
   char*     class_filename,
   double    positive_constraint,
   double    negative_constraint,
   double    convergence_threshold,
   long      seed,
   long int  iterations,
   double    objective,
   BOOLEAN_T print_time,
   FILE*     outfile)
{
  /* Print the name of the program. */
  if (COMPUTE_WEIGHTS) {
    fprintf(outfile, "# Generated by compute-weights\n");
  } else if (KERNEL_PCA) {
    fprintf(outfile, "# Generated by kernel-pca\n");
  } else if (KFD) {
    fprintf(outfile, "# Generated by kfd\n");
  }

  // Print release, web site and citation info.
  fprintf(outfile, "# Gist, version %s\n", RELEASE);
  fprintf(outfile, "# For more information, go to http://svm.sdsc.edu\n# \n");
  // FIXME: Add cite.
  // fprintf(outfile, "# If you use this software in your research, please cite:\n");

  /* Print the required inputs. */
  fprintf(outfile, "# train_file=%s\n", kernel->train_filename);
  if (class_filename == NULL) {
    fprintf(outfile, "# class_file=NULL\n");
  } else {
    fprintf(outfile, "# class_file=%s\n", class_filename);
  }

  /* Print the properties of the kernel. */
  fprintf(outfile, "# matrix_from_file=%s", 
	  boolean_to_string(kernel->matrix_from_file));
  fprintf(outfile, " zero_mean=%s variance_one=%s", 
	  boolean_to_string(kernel->zero_mean_row),
	  boolean_to_string(kernel->variance_one));
  fprintf(outfile, "  normalize=%s\n",
	  boolean_to_string(kernel->normalize));
  fprintf(outfile, "# constant=%g", kernel->constant);
  fprintf(outfile, " coefficient=%g", kernel->coefficient);
  fprintf(outfile, " power=%g\n", kernel->power);
  fprintf(outfile, "# radial=%s", boolean_to_string(kernel->radial));
  fprintf(outfile, " width_factor=%g", kernel->width_factor);
  fprintf(outfile, " two_squared_width=%g", kernel->two_squared_width);
  fprintf(outfile, " add_diag=%g\n", kernel->add_diag);

  if (COMPUTE_WEIGHTS) {

    // Print feature selection parameters.
    fprintf(outfile, "# feature_select=%s", 
	    convert_enum_type(kernel->feature_select, FSELECT_STRS, 
			      NUM_FSELECT_T));
    fprintf(outfile, " thresh_type=%s",
	    convert_enum_type(kernel->thresh_type, THRESH_STRS, NUM_THRESH_T));
    fprintf(outfile, " fthreshold=%g\n", kernel->fthreshold);

    // Print the properties of the training.
    fprintf(outfile, "# positive_constraint=%g negative_constraint=%g",
	    positive_constraint, negative_constraint);
    fprintf(outfile, " constrain_weights=%s\n", 
	    boolean_to_string((positive_constraint != 0.0) ||
			      (negative_constraint != 0.0)));
    fprintf(outfile, "# positive_diagonal=%g negative_diagonal=%g\n",
	    kernel->positive_diagonal, kernel->negative_diagonal);
    fprintf(outfile, "# convergence_threshold=%g seed=%ld\n",
	    convergence_threshold, seed);
    fprintf(outfile, "# objective=%g iterations=%ld\n", objective, iterations);
  }
  if (KFD) {
    fprintf(outfile, "# regularizer=%g\n", kernel->regularizer);
  } else {
    kernel->regularizer = 0.0;
  }

  if (print_time) {
    fprintf(outfile, "# time=%.8g s", myclock() / 1.0E6);
    fprintf(outfile, " host=%s", hostname());
    fprintf(outfile, " date=%s\n", date_and_time());
  }
}

#ifdef MAIN
#include "cmdline.h"

#define CONVERGENCE_THRESHOLD 0.000001  // Minimum delta in objective function.
#define MAXITER 10000                   // Maximum number of iterations.
#define EIGENVALUE_THRESHOLD 0          // Default threshold for eignevalues.
#define MINIMUM_PRECISION 4             // At least this many digits in output.

VERBOSE_T verbosity = NORMAL_VERBOSE;

int main(int argc, char *argv[]) {

  char*     class_filename = NULL;        // Classifications of training data.
  FILE*     class_file = NULL;
  int       class_number = 0;             // In a multiclass file, which class to use. First one is 1

  // Compute weights options.
  double    hold_out = 0.0;               // Percent to do hold-one-out.
  char*     score_filename = NULL;        // Feature selection score output.
  char*     initial_filename = NULL;      // Initial weights file.
  double    positive_constraint = 0.0;    // Maximum weight for positives.
  double    negative_constraint = 0.0;    // Maximum weight for negatives.
                                          // Delta objective at which to stop.
  double    convergence_threshold = CONVERGENCE_THRESHOLD; 
  long int  maxiter = MAXITER;            // Maximum number of iterations.

  // Kernel PCA options.
  int       num_eigens = 0;               // Number of eigenvectors to produce.
  double    eigen_threshold = EIGENVALUE_THRESHOLD;  // Minimum eigenvalue.
  char*     eigenvalue_filename = NULL;   // File to print eigenvalues.

  // Generic options.
  BOOLEAN_T format_line = FALSE;          // Use format line in RDB files?
  BOOLEAN_T kernel_output = FALSE;        // Output kernel matrix and stop?
  long      seed = time(0);               // Seed for random number generator.
  BOOLEAN_T print_time = TRUE;            // Include timing data in output?

  // Data structures.
  KERNEL_T*      kernel = allocate_kernel(TRUE); // Kernel function.
  CLASS_ARRAY_T* classes = NULL;          // The training set classifications.
  ARRAY_T*       classes_array = NULL;    // Same as above, but raw array.
  ARRAY_T*       score_array = NULL;      // Feature quality scores.
  ARRAY_T*       weights = NULL;          // The weights.
  ARRAY_T*       classifications = NULL;  // Predicted classifications.
  ARRAY_T*       discriminants = NULL;    // Training set discriminants.
  MATRIX_T*      output_matrix = NULL;    // Weights OR principal components.
  RDB_MATRIX_T*  rdb_output_matrix = NULL;// RDB version of previous.

  // Local variables.
  long int       iterations = 0;    // Total number of training iterations.
  double         objective = 0.0;   // Final value of objective function.
  int            output_precision;  // Number of digits in output.
  verbosity = NORMAL_VERBOSE;

  // Be sure to zero-mean the columns for kernel PCA.
  if (KERNEL_PCA) {
    kernel->zero_mean_col = TRUE;
  }

  // Parse the command line.
  DO_STANDARD_COMMAND_LINE
    (2,
     DATA_OPTN(1,
	       train, <filename> (required), 
	       kernel->train_filename = _OPTION_);
     DATA_OPTN(COMPUTE_WEIGHTS + KFD,
	       class, <filename> (required), class_filename = _OPTION_);
     DATA_OPTN(COMPUTE_WEIGHTS, useclassnumber, <value>, 
	       class_number = atoi(_OPTION_));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       initial, <filename>, initial_filename = _OPTION_);
     FLAG_OPTN(1,
	       matrix, kernel->matrix_from_file = TRUE);
     DATA_OPTN(COMPUTE_WEIGHTS,
	       holdout, <percent>, hold_out = atof(_OPTION_) / 100.0);
     FLAG_OPTN(1,
	       zeromeanrow, kernel->zero_mean_row = TRUE);
     FLAG_OPTN(1,
	       varone, kernel->variance_one = TRUE);
     DATA_OPTN(COMPUTE_WEIGHTS,
	       //	       fselect, fisher|ttest|welch|mannwhitney|sam|samf|tnom,
	       fselect, fisher|ttest|welch|mannwhitney|sam|tnom, // hide the samf option for now.
	       kernel->feature_select = convert_enum_type_str(_OPTION_,
							      INVALID_FSELECT,
							      FSELECT_STRS,
							      NUM_FSELECT_T));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       fthreshtype, percent|number|value,
	       kernel->thresh_type = convert_enum_type_str(_OPTION_,
							   INVALID_THRESH,
							   THRESH_STRS,
							   NUM_THRESH_T));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       fthreshold, <value>, kernel->fthreshold = atof(_OPTION_));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       fscores, <file>, score_filename = _OPTION_);
     FLAG_OPTN(1,
	       nonormalize, kernel->normalize = FALSE);
     DATA_OPTN(1,
	       constant, <value> (default=1), 
	       kernel->constant = atof(_OPTION_));
     DATA_OPTN(1,
	       coefficient, <value> (default=1), 
	       kernel->coefficient = atof(_OPTION_));
     DATA_OPTN(1,
	       power, <value> (default=1), kernel->power = atof(_OPTION_));
     FLAG_OPTN(1,
	       radial, kernel->radial = TRUE);
     DATA_OPTN(1,
	       widthfactor, <value> (default = 1),
     	       kernel->width_factor = atof(_OPTION_));
     DATA_OPTN(1,
	       adddiag, <value> (default=0), 
	       kernel->add_diag = atof(_OPTION_));
     FLAG_OPTN(KERNEL_PCA,
	       nocenter, kernel->zero_mean_col = FALSE);
     DATA_OPTN(KERNEL_PCA,
	       numeigens, <value> (default = all),
	       num_eigens = atoi(_OPTION_));
     DATA_OPTN(KERNEL_PCA,
	       eigenthresh, <value> (default = 0.000001),
	       eigen_threshold = atof(_OPTION_));
     DATA_OPTN(KERNEL_PCA,
	       eigenvalues, <file>, eigenvalue_filename = _OPTION_);
     DATA_OPTN(KFD,
	       regularizer, <value> (default = 1),
	       kernel->regularizer = atof(_OPTION_));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       posconstraint, <value> (default = 1), 
	       positive_constraint = atof(_OPTION_));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       negconstraint, <value> (default = 1), 
	       negative_constraint = atof(_OPTION_));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       diagfactor, <value> (default = 0.1),
	       kernel->diagonal_factor = atof(_OPTION_));
     FLAG_OPTN(1,
	       rdb, format_line = TRUE);
     FLAG_OPTN(1,
	       kernelout, kernel_output = TRUE);
     DATA_OPTN(COMPUTE_WEIGHTS,
	       threshold, <value> (default = 0.00001),
	       convergence_threshold = atof(_OPTION_));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       maxiter, <value> (default = 10000),
	       maxiter = atoi(_OPTION_));
     DATA_OPTN(COMPUTE_WEIGHTS,
	       seed, <value> (default from clock),
	       seed = atoi(_OPTION_));
     FLAG_OPTN(1,
	       notime, print_time = FALSE);
     DATA_OPTN(1,
	       verbose, 1|2|3|4|5 (default=2),
	       verbosity = (VERBOSE_T)atoi(_OPTION_));
     );

  // Initialize the random number generator.
  my_srand(seed);

  // Initialize the timer.
  myclock();

  // Print the process ID.
  if (verbosity >= NORMAL_VERBOSE) {
    fprintf(stderr, "process_id=%d\n", (int)(getpid()));
  }

  // The output precision depends on the convergence threshold.
  output_precision = (int)(abs(log10(convergence_threshold))) - 2;
  if (output_precision < MINIMUM_PRECISION) {
    output_precision = MINIMUM_PRECISION;
  }

  // Set the type of feature selection.
  set_default_feature_selection(kernel);

  // Make sure we got the required files.
  if (kernel->train_filename == NULL) {
    die("No training file given.");
  }
  if (!KERNEL_PCA) {
    if (class_filename == NULL) {
      die("No classification file given.");
    }
    if (class_number < 0) { // this allows -classnumber 0, which gives the same behavior as -classnumber 1 and the default behavior.
      die("-classnumber must be invoked with a value of 1 or more.");
    }
  }

  // Make sure we don't have conflicting options.
  if (kernel->zero_mean_row && kernel->matrix_from_file) {
    die("Cannot use zero_mean option with kernel matrix input.");
  }
  if (kernel->variance_one && kernel->matrix_from_file) {
    die("Cannot use variance_one option with kernel matrix input.");
  }
  if ((kernel->feature_select != NONE_FSELECT) && kernel->matrix_from_file) {
    die("Cannot use feature selection with kernel matrix input.");
  }
  if ((kernel->feature_select == NONE_FSELECT) && score_filename) {
    die("Cannot produce a feature score file without feature selection.");
  }
  if (((positive_constraint != 0.0) || (negative_constraint != 0.0)) &&
      (kernel->diagonal_factor != 0.0)) {
    fprintf(stderr, "Warning: Using 1-norm and 2-norm soft margins simultaneously.\n");
    fprintf(stderr, "1-norm box constraint: (%g,%g)\n", positive_constraint,
	    negative_constraint);
    fprintf(stderr, "2-norm diagonal factor: %g\n", kernel->diagonal_factor);
  }

  // Read the classifications.
  if (!KERNEL_PCA) {
    if (open_file(class_filename, "r", TRUE, "class", "classes",
		  &class_file) == 0)
      exit(1);
    classes = read_classifications(format_line, class_file, class_number);
  } else {
    // Can't have a diagonal factor with kernel PCA (no classes).
    kernel->diagonal_factor = 0.0;
  }

  // Read the data files.
  read_data_files(format_line, kernel);

  // Compute the base kernel matrix.
  compute_base_kernel(kernel, classes, &score_array);
  
  // Turbo-ize the kernel matrix.
  transform_kernel_matrix(TRUE,
			  kernel, 
			  classes, 
			  kernel->self_train,
			  kernel->self_test,
			  kernel->kernel_matrix);

  // Stop here if requested.
  if (kernel_output) {
    print_rdb_matrix(kernel->kernel_matrix_rdb, format_line, 8, 6, stdout);
    exit(0);
  }

  if (COMPUTE_WEIGHTS) {

    // Allocate and initialize the weights array.
    weights = initialize_weights(format_line, 
				 get_num_rows(kernel->kernel_matrix),
				 initial_filename);

    // Optimize the weights.
    iterations = double_optimize_weights(convergence_threshold, 
					 maxiter,
					 positive_constraint,
					 negative_constraint,
					 classes,
					 kernel->kernel_matrix,
					 weights);

    // Compute the final value of the objective function.
    objective = compute_objective(weights,
				  classes, 
				  kernel->kernel_matrix);

    // Encode the classifications as the signs of the weights.
    sign_weights(classes, weights);

    // Create the output matrix. 
    if (hold_out > 0.0) {
      output_matrix = allocate_matrix(get_array_length(weights), 
				      NUM_COLUMNS);
    } else {
      output_matrix = allocate_matrix(get_array_length(weights), 
				      NUM_COLUMNS - 2);
    }

    // Store the classes in the matrix.
    classes_array = get_class_array(classes);
    set_matrix_column(classes_array, TRUE_CLASS_COLUMN, output_matrix);

    // Store the weights in the matrix.
    set_matrix_column(weights, WEIGHT_COLUMN, output_matrix);

    // Classify the training set.
    classify_list(get_num_rows(kernel->kernel_matrix),
		  weights,
		  kernel->kernel_matrix,
		  &classifications,
		  &discriminants);
    set_matrix_column(classifications, TRAIN_CLASS_COLUMN, output_matrix);
    set_matrix_column(discriminants, TRAIN_DISC_COLUMN, output_matrix);

    // If requested, do hold-one-out cross-validation.
    if (hold_out > 0.0) {
      free_array(classifications);
      free_array(discriminants);
      hold_out_classify_list(hold_out,
			     convergence_threshold, 
			     maxiter,
			     positive_constraint,
			     negative_constraint,
			     format_line,
			     score_filename,
			     kernel,
			     classes,
			     weights, 
			     score_array,
			     &classifications,
			     &discriminants);

      set_matrix_column(classifications, HOLD_OUT_CLASS_COLUMN, output_matrix);
      set_matrix_column(discriminants, HOLD_OUT_DISC_COLUMN, output_matrix);
    }
  }


  if (KERNEL_PCA) {
    // Compute the normalized eigenvectors.
    if (verbosity >= NORMAL_VERBOSE) {
      fprintf(stderr, "Computing eigenvectors of the kernel matrix.\n");
    }
    output_matrix = find_normalized_eigenvectors(num_eigens, 
						 eigen_threshold,
						 eigenvalue_filename,
						 kernel->kernel_matrix);
  }

  if (KFD) {
    // Find the kernel fisher discriminant.
    output_matrix = find_kfd(kernel->regularizer,
			     classes,
			     kernel->kernel_matrix);
  }

  // Create an RDB version of the output matrix.
  rdb_output_matrix = allocate_rdb_matrix(get_num_rows(output_matrix), 
					  get_num_cols(output_matrix),
					  output_matrix);
  set_row_names(get_row_names(kernel->kernel_matrix_rdb), rdb_output_matrix);
  set_corner_string(get_corner_string(kernel->kernel_matrix_rdb),
		    rdb_output_matrix);

  if (COMPUTE_WEIGHTS) {
    add_col_name("class", rdb_output_matrix);
    add_col_name("weight", rdb_output_matrix);
    add_col_name("train_classification", rdb_output_matrix);
    add_col_name("train_discriminant", rdb_output_matrix);
    if (hold_out > 0.0) {
      add_col_name("hold_out_classification", rdb_output_matrix);
      add_col_name("hold_out_discriminant", rdb_output_matrix);
    }
  }

  if (KERNEL_PCA) {
    add_eigen_col_names(rdb_output_matrix);
  }
  if (KFD) {
    add_col_name("kfd", rdb_output_matrix);
  }

  // Write all the parameters to the header.
  print_weight_header(kernel,
		      class_filename,
		      positive_constraint,
		      negative_constraint,
		      convergence_threshold,
		      seed,
		      iterations,
		      objective,
		      print_time,
		      stdout);

  // Print the weights OR eigenvectors.
  print_rdb_matrix(rdb_output_matrix, format_line, output_precision + 2, 
		   output_precision, stdout);
#ifdef DEBUG
  // Free dynamic memory.
  free_kernel(kernel);
  free_class_array(classes);
  free_array(classes_array);
  //free_array(score_array); Don't want to free multiple times.
  free_array(weights);
  free_array(classifications);
  free_array(discriminants);
  free_rdb_matrix(rdb_output_matrix);
#endif

  return(0);
}
#endif


/*
 * Local Variables:
 * mode: c
 * c-basic-offset: 2
 * End:
 */
